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FOREWORD 


This  final  report  was  prepared  by  the  Aerodynamics  Section  of  the  Grumman 
Aerospace  Corporation,  Bethpage,  New  York  for  the  Flight  Mechanics  Division, 

Air  Force  Flight  Itynamics  Laboratory,  Wright-Patterson  Air  Force  Base,  Ohio. 

The  work  was  performed  under  Contract  No.  F33615-75-C-3073,  which  was  initiated 
under  Project  No.  1476,  "Advanced  Wing-Body  Aerodynamic  Analysis  and  Design." 

Mr.  J.  Kenneth  Johnson  (FXM)  was  the  Project  Monitor  of  this  contract. 

The  report  consists  of  three  volumes.  Volume  I,  entitled  "Description  of 
Analysis  Methods  and  Applications,"  describes  the  methods  used  to  predict  sur- 
face pressure  distributions  and  aerodynamic  forces  on  three-dimensional  wing- 
body  combinations  at  transonic  speeds,  including  viscous  effects.  Volume  I 
also  contains  an  extensive  set  of  comparisons  between  numerical  predictions 
and  experimental  results.  Detailed  instructions  required  to  use  the  program 
are  provided  in  Volume  II,  "User's  Manual  and  Code  Description."  Volume  III 
was  written  at  Sybucon,  Inc.,  Atlanta,  Georgia  and  contains  a complete  descrip- 
tion of  the  theory  and  program  that  computes  the  full  three-dimensional  boundary 
layer  over  the  wing.  This  work  was  performed  by  Sybucon  under  subcontract  to 
Grumman  Aerospace.  Although  this  program  operates  independently  of  the  program 
described  in  Volume  II,  the  input  data  set  required  for  the  full  three-dimensional 
boundary  layer  computation  is  generated  by  the  code  documented  in  Volume  II. 

Mr.  F.  Berger  was  the  Program  Manager;  Dr.  W.  Mason  and  Mr.  D.  Mackenzie 
served  as  Project  Engineers.  The  work  was  performed  in  close  cooperation  with 
the  co-authors  from  the  NASA  Ames  Research  Center,  Dr.  W.  F.  Ballhaus  and 
Ms.  J.  Frick.  Additional  contributors  to  the  project  included  G.  Simpers, 

A.  Vachris,  D.  Raila,  P.  Aidala,  M.  Sturm  and  A.  Bunnell  of  Grumman,  and 
Drs.  F.  R.  Bailey  and  T.  Holst  of  NASA  Ames.  Moreover,  contributions  have  been 
made  by  A.  Chen  of  Boeing,  Drs.  R.  Melnik,  B.  Grossman  and  G.  Volpe  of  the 
Grumman  Research  Department  and  Grumman  Consultants  Prof.  A.  Jameson,  Prof.  J. 
Werner  and  Dr.  E.  Murman.  As  noted  above,  the  three-dimensional  boundary  layer 
program  was  written  by  Dr.  J.  Nash  and  Dr,  R.  Scruggs  of  Sybucon,  Inc. 
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SECTION  I 


INTRODUCTION 

1.  BACKGROUND 

Hie  method  described  herein  was  developed  with  the  objective  of 
producing  a computationally  efficient  scheme  for  calculating  three- 
dimensional  boundary  layers:  laminar  or  turbulent,  over  a finite  swept 
wing  of  arbitrary  planform  and  thickness,  in  compressible  steady  flow. 

The  method  is  intended  to  be  used  both  for  routine  boundary- layer  calcu- 
lations in  specified  pressure  fields,  and  for  interactive  calculations 
in  which  the  pressure  field  will  be  generated  by  a parallel  potential 
flow  calculation  via  suitable  iterative  coupling. 

The  approach  adopted  has  been  to  draw  on  the  extensive  experience 
gained  with  the  earlier  turbulent  methodology  (References  [1  through  81), 
to  retain  the  turbulence-model  formulation  essentially  intact,  and  to 
concentrate  on  improving  the  numerical  aspects  of  the  method  and  its 
data-handling  capabilities.  The  need  for  radical  improvements  in  this 
latter  area  became  evident  from  the  large  amount  of  hand  work  which  had 
to  be  done  in  performing  calculations  even  for  the  simple  planform  geom- 
etry considered  in  References  [6,7].  It  was  subsequently  decided  to 
upgrade  the  basic  fluid-dynamic  package  by  including  the  thermal  energy 
equation,  in  differential  form,  among  the  governing  equations,  so  as  to 
provide  greater  generality  for  the  method.  It  was  also  decided  to  avoid 
the  limitations,  imposed  by  matching  the  solution  to  a separate  inner- 
layer  calculation  near  the  wall,  by  continuing  the  numerical  solution 
through  the  viscous  sublayer.  Lastly,  the  provision  for  performing  purely 
laminar  calculations,  near  the  leading  edge  or  elsewhere,  was  to  be  made. 


N 


As  far  as  Che  numerical  scheme  was  concerned,  It  was  thought  essen- 
tial to  replace  the  explicit  formulation,  used  In  the  earlier  methods, 
by  an  Implicit  one.  Only  by  doing  so  would  it  be  possible  for  the  user  to 
gain  control  over  the  chordwlse  step  length  and  escape  from  the  excessive 
computer  run-times  which  had  been  typical  at  high  Reynolds  numbers.  On 
the  other  hand,  the  adoption  of  an  implicit  numerical  scheme  posed  a num- 
ber of  problems  which  had  not  arisen  with  the  older  explicit  method. 

First,  the  long  chordwlse  step  lengths,  maae  possible  by  the  implicit 
scheme,  and  the  need  to  scan  the  integration  domain  laterally,  presented 

problems  in  proper  treatment  of  the  zones  of  dependence.  These  problems 

were  resolved  by  using  an  alternating  direction  technique,  which  offered 
other  advantages  too.  Second,  the  implicit  formulation  of  the  particular 
set  of  governing  equations,  which  Included  the  empirically-modified 
turbulent  kinetic-energy  equation,  posed  stability  problems  associated 
with  changing  character  of  the  equations:  the  equations  become  hyperbolic 
near  the  wall  whereas  (since  molecular  viscosity  is  included)  they  are 
parabolic  in  the  outer  part  of  the  flow.  The  stability  problems  were 
resolved  by  representing  the  turbulent  dissipation  as  a sink  term  in  the 
matrix  form  of  the  governing  equations. 

2.  OUTLINE  OF  THE  METHOD 

The  wing  surface  is  segmented  into  quasi-quadrilateral  panels 
(Figure  1)  and  a locally  orthogonal  set  of  curvilinear  coordinates  is 
constructed  following  a conical  development  of  the  local  region.  The 
development  is  required  to  satisfy  certain  geometrical  constraints. 
Streamwise  cuts  are  made  over  the  planform  at  arbitrary  spanwise  loca- 
tions and  the  number  of  surface  points  on  each  cut  is  the  same.  The 

2 


J 


calculation  then  proceeds  from  the  planform  leading  edge;  proper  account 
being  taken  as  to  whether  the  leading  edge  is  turbulent  or  laminar.  The 
governing  equations  are  f inite-dif ferenced  with  respect  to  the  local 
surface  normal  coordinate  y,  and  the  locally  transverse  coordinate  z,  so 
that  in  the  local  y-z  plane  'he  difference  equations  are  implicit.  The 
solution  is  then  obtained  at  each  x (streamwise)  station  in  succession. 
This  chordwise  forward-marching  is  applied  with  arbitrary  step  length 
proceeding  from  leading  to  trailing  edge  or  until  separation  is  encoun- 
tered at  some  spanwise  location.  Provision  is  made  to  step  past  locally 
separating  regions  if  desired.  Transition  from  laminar  to  turbulent  flow 
may  be  triggered  at  the  user's  discretion  at  any  chordwise  position, 
except  when  the  leading  edge  flow  is  already  turbulent. 
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U u,  V + v,  VI  + w,  respectively,  where  U,  V,  W are  defined  as  time- 
averaged  means,  and  u,  v,  w are  the  Instantaneous  fluctuations  about 


E 


thoae  means.  Similarly,  the  static  temperature,  pressure,  and  density 
are  expressed  as  T + ♦,  P + p' , P + p* , respectively. 

The  equations  of  motion:  consisting  of  the  two  mean-flow  momen- 
tum equations,  the  mean  thermal-energy  equation,  and  the  continuity 
equation,  can  accordingly  be  written  in  the  form 


£U  3U  + _1_ 
hj  3x  hj 


(pV  + p'v) 


3U 

3y 


£W 

b. 


3U 

3z 


+ pW(KnU  - K31W)  + 


_1_ 

hi  3x 


13.-—.  1*  ,j l 31). 

+ ^ ~ (puv)  - - — (^ 
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en  AW  + A.  (pV  + rrv)  aw  + pw  aw 

hj  3x  i>2  3y  KT  3z 
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+ pu(K21  + K31)  + pW(K13  + K23)  - 0 , 


(2.9) 
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where  c is  Che  race  of  Curbulenc  dissipacion  inCo  heaC. 

In  Equations  (2.6,  7,  8),  p,  k and  are  Che  molecular  viscosiCy, 
Che  Chermal  conductivity,  and  Che  specific  heat  at  constant  pressure, 
respectively,  and  they  can  be  related  to  one  another  via  the  Prandtl  num- 
ber: Pr  ■ pc  /k,  which  is  taken  to  be  a constant  (-  0.72);  c is  also 
P P 

2 2 

taken  to  be  a constant  (»  6006  ft  .sec  °R) , and  hence  p and  k are  pro- 
portional to  each  other.  The  molecular  viscosity  is  taken  to  be  a func- 
tion of  local  static  temperature  via  the  Sutherland  formula: 


3/2 

p - 2.270  t -7  X9876  x 10~8  lb  sec/ft2  • 


(2.10) 


The  fluid  is  assumed  to  be  thermally  perfect,  whereupon  the  equation  of 
state  holds,  in  the  form 


p - pRT  (2.11) 

2 2 

where  the  gas  constant:  R - 1716  ft  /sec  °R.  These  values  of  c and 

r 

R are  consistent  with  y * 1.4,  where  y is  the  ratio  of  the  specific  heats. 
3.  CLOSURE  RELATIONSHIPS 

Closure  of  the  system  of  equations  is  effected  by  means  of  the 
turbulent  kinetic-energy  equation,  which  is  written  in  the  form 


(2.12) 


i 


sS 


where 


2 2 2 2 
q - u + v + w 


(2.13) 


together  with  empirical  relationships  to  model  the  diffusion  and  dissi- 
pation terms,  and  to  model  the  dependence  of  the  turbulent  shear 

~2 

stresses  and  turbulent  heat  flux  on  q . 

Specifically,  following  Bradshaw  et  al.  [10,  11]  and  Nash  [1],  the 
diffusion  and  dissipation  terms  are  modeled,  respectively  by 


£_v  . 5_v 
2 


max  2 

Q~  a2  q 
e 


(a2)372 


(2.14) 


(2.15) 
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where  q is  the  maximum  value  of  q in  the  outer  three-fourths  of  the 

max 

boundary- layer,  and  where  the  diffusion  function:  a^,  and  the  dissipa- 
tion length:  L^,  are  assumed  to  be  universal  functions  of  physical  dis- 
tance through  the  boundary  layer.  The  particular  forms  for  a ^ and  L^, 
used  here,  are  the  same  as  those  found  appropriate  in  earlier  studies 
[3,4],  and  can  be  represented  by 


a2  - 1.125  n2  - 0.375  n4 


7.195 


2 7 

1 + 4n  + 5n 


(2.16) 


(2.17) 


These  two  functions  are  illustrated  in  Figure  3.  In  Equations  (2.16,17) 
q is  the  dimensionless  physical  distance  through  the  boundary  layer: 


1 a 

n " 6 J h2dy  * 


(2.18) 
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and  6 is  Che  boundary- layer  thickness,  defined  as  Che  distance  from  the 
surface  where 


«Ue- 


ur  + (w 


- W)2} 


1/2 


0.005  Q 


(2.19) 


The  turbulent  shear-stress  components  are  related  to  q by  empiri- 
cal functions  of  the  form 


puv 


(2.20) 


pvw 


(2.21) 


Several  earlier  models  are  contained  as  special  cases  in  Equations 

(2.20,21):  If  f(^)  = constant,  the  Townsend  structural-equilibrium  model 

[12]  is  recovered;  if  f(i ]))  * i|/,  the  isotropic  eddy  viscosity  model  is 

recovered  (although  not  a prescribed  eddy  viscosity)  finally,  if 
2 

f(i /')*</'  , the  mixing  length  model  is  recovered.  It  is  the  intention  to 
explore  the  implications  of  and  experimental  support  for,  alternative 
forms  for  the  function  f at  a later  date;  meanwhile,  a cursory  study  of 
some  available  experimental  data  (Figure  4)  indicates  the  suitability  of 
the  form 

f(<J»)  - 0.0225  <|i,  (2.22) 

and  this  form  has  been  provisionally  adopted. 

~2 

The  corresponding  relationship  between  turbulent  heat  flux  and  q 
is  written  as 
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(2.23) 


where  g is  another  empirical  function.  To  date  the  various  possible 
forms  for  g have  not  been  explored,  and  the  present  work  is  based  on 
the  provisional  form 


g(^T) 


0.0225  ik 


(2.24) 


in  which  Prfc  is  an  assumed  constant  value  (»  0.9)  of  the  turbulent 
Prandtl  number.  An  obvious  refinement  could  be  to  allow  Prt  to  vary 
with  position  across  the  boundary  layer. 

4.  DIMENSIONLESS  VARIABLES 

Characteristic  scales  of  length,  velocity,  temperature,  and  density: 

L , Q , T , p , respectively,  are  introduced  for  the  purposes  of  reducing 
o o o o 

the  variables  to  dimensionless  form.  Two  of  these,  in  turn,  define  a 

Mach  number,  M , where 
o 


M Try 

° (pRTo)1/2 


(2.25) 


the  particular  choice  of  characteristic  scales  is  of  no  special  signifi- 
cance; here,  Q is  identified  with  the  free-stream  velocity,  and  T and 
o o 

p are  associated  with  stagnation  conditions, 
o 

The  dimensionless  variables  subsequently  become 


x - y - z 

L ’ y L * L 
o o o 


(2.26) 
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2 — - — 

2 q — uv  — vw 

q ’ 72  • uv  " 7T  • w " 72 

Qo  Qo  *0 


Q T * 
xo  o 


(2.27) 


(2.28) 


(2.29) 


(2.30) 


(2.31) 


In  what  follows,  the  - sign  will  be  dropped  for  clarity,  but  it  will  be 
understood  that  the  variables  retain  their  dimensionless  status  as 
defined  here. 

5*  EQUATIONS  IN  MATRIX-VECTOR  FORM 

The  two  momentum  equations  (Equations  2.6,7),  the  thermal-energy 

equation  (Equation  2.8),  and  the  empirically  modified  turbulent  kinetic- 

energy  equation  (Equation  2.12)  are  selected  as  the  principle  governing 

equations,  leaving  the  continuity  equation  to  be  solved  retrospectively. 


These  four  principle  governing  equations  can  conveniently  be  written  in 
matrix  vector  form,  as  follows: 


A F + A — + A — ^ + A 
V A2  3x  A3  A4  3y 


32F 

* *5  72  + *6  ■ ° • 
3y 


(2.32) 


where  F is  a four-dimensional  vector  of  the  principle  dependent  varla- 


FT  - {U,  W,  T,  q2}  . 


(2.33) 


In  Equation  (2.32)  the  A3  through  A^  are  4x4  matrices,  and  Afc  is  a 
4-dimensional  vector. 

and  A^  are  simply  products  of  a scalar  with  the  identity 
matrix,  I,  of  order  4: 


a ■ t a * I • 
2 hx  ’ 3 h3  * 


(2.34) 


A,.  is  a diagonal  matrix,  and  A ^ is  full.  A^  could  contain  off-diagonal 
terms,  depending  on  the  particular  formulation;  here  it  contains  only 
one  term,  as  indicated  below. 


Writing 


A1  " ^(Di.j1 


h2A4  ‘ (a(4)i.J1 


h2A5  " la(5)i, j1 


A6  “ (*(6)i}  ’ 


the  non-zero  content  of  these  matrices  is  as  follows: 
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a(4)4,l  " 2p  q fu  ’ 


a(4)4,2  " 2p  q fw  * 


2pW2 


a(4)4.4  " pV  4 0 


*(5)1,1  “ '(l‘  + "u)  * 


*(5)2,2 
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(2.41) 


(2.42) 


(2.43) 


(2.44) 


(2.45) 


(2.46) 


(2.47) 


(2.48) 


(2.49) 


(2.50) 


(2.51) 


(2.52) 


(2.53) 


. 1 Ay. 

U uh2  3y  * 


(2.63) 


1 3k 
Ck  kh2  3y 


(2.64) 


In  the  above  equations,  f',  g'  denote  differentiation  with  respect  to 
their  arguments;  the  subscripts:  u,  w are  attached  to  f,  f'  to  distin- 
guish between  the  arguments  involving  3U/3y  and  3W/3y,  respectively. 

6.  BOUNDARY  CONDITIONS 

The  boundary  conditions  on  the  surface  over  which  the  boundary 
layer  is  developing,  which  is  assumed  to  be  a smooth,  impervious,  ther- 
mally insulated  wall,  are: 


U - V - W - 0 (2.65) 

3T/3y  - 0 (2.66) 

q2  - 0 (2.67) 

At  the  outer  edge  of  the  boundary  layer  the  appropriate  conditions  are: 


. i”  . U - o 

3y  3y  3y 


(2.68) 


(2.69) 


Initial  conditions,  at  the  leading  edge,  are  assumed  to  conform  to  the 
similarity  relationships  pertaining  to  asymptotic,  attachment- line  flow 
at  the  local  leading-edge  sweep  angle.  Under  such  conditions  the  bound- 
ary layer  quantities  may  be  expressed  in  the  form 
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I* 


F - F(y,  C*,  Me)  , (2.70) 

where  F is  the  vector  dependent  variable  (Equation  (2.33),  M is  the 

e 

local  external  Mach  number,  and  C*  is  the  similarity  parameter 

C*  - p W2  /(Vi  ai:  /Jx  ) , (2.71) 

we  we  n 
P » 

in  which  W is  the  external  veiociiv  component  parallel  to  the  attach- 

CP 

ment  line  U is  the  component  normal  to  it,  and  x is  the  normal  dis- 
e n 

n 

tance  from  the  attachment  line,  measured  along  the  wing  contour.  The 
flow  in  the  boundary  layer,  on  the  attachment  line,  is  considered  to 
be  laminar  (at  least  for  » 0)  if  C*  < 10^  (see  Reference  [13]),  and 
turbulent  if  C*  > 10^. 

In  the  present  work,  the  functional  form,  corresponding  to  Equa- 
tion (2.70),  is  initially  taken  to  be  the  same  as  its  incompressible 
equivalent,  and  subsequently  is  improved  upon  via  an  iteration  involving 
conditions  at  the  first  two  chordwise  stations.  The  incompressible 
form  for  laminar  flow  is  tabulated  in  the  text  books  (see,  e.g..  Refer- 
ence [14]),  and  experimental  results  for  incompressible  flow  are  given 
in  References  [13]). 

Boundary  conditions  must  also  be  specified  along  the  sides  of 
the  Integration  domain:  i.e.  for  the  extreme  values  of  z,  which  corre- 
spond to  the  root  and  tip  of  the  wing.  The  true  conditions  along  these 
boundaries  are,  of  course,  unspecified  and  the  most  that  can  be  done  is 
to  provide  the  user  with  options,  representing  stereotyped  boundary  con- 
ditions, permitting  sensitivity  studies  to  be  carried  out  for  the 
particular  configuration  and  pressure  distribution.  Candidate  stereotyped 
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root  and  tip  conditions  considered  here  are: 

(a)  two-dimensional  flow 

(b)  plane-of-symmetry  flow 

(c)  infinite-yawed-wing  flow. 

In  each  case,  different  assumptions  are  made  concerning  the 
spanwise  derivatives  on  the  boundaries,  and  these  assumptions  take  the 
place  of  detailed  information  about  the  flow  quantities  being  convected 
across  those  boundaries. 

7.  DISPLACEMENT  THICKNESS  AND  WALL  SHEAR  STRESS 

The  displacement  thickness:  5*,  in  a three-dimensional  flow 
(Reference  [4,15)),  can  be  determined  by  matching  the  normal  component 
of  mass  flux,  in  the  real  flow,  to  that  which  would  exist  in  a hypotheti- 
cal inviscid  flow  about  a displacement  surface  whose  distance  from  the 
wall  is  5*.  In  this  equivalent  inviscid  flow: 


PeUe  35*  . peWe  36* 


h^  3x 


h^  3z 


+ / 3?  (peVe)dy  * 

yD 


(2.72) 


where  3(p  V )/3y  is  given  by  the  continuity  equation  (Equation  2.9): 
e e 


r2  k ‘‘’.V  ' \ W + ry  li  ‘"eV 


+ »eUe(K21  + K31}  + peVK13  + K23>  ‘ (2'73) 


and  y^  is  the  value  of  y corresponding  to  6*.  Upon  substitution  into 
Equation  (2.72),  we  have 
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,V  + / (peUeKn  + Pe»eK23)h2dy 


PeUe  36*  + PeWe  36* 
h,  3x  h_  3z 


— r~-  ~ — (p  U ) + r~“  ~ — (p  W ) 
hj  3x  e e 3z  e e 


+ peUeK31  + peWeK32  / h2dy 

yD 


(2.74) 


in  which  the  integral  on  the  right-hand  side  represents  the  physical 
distance  between  and  y.  The  left-hand  side  is  set  equal  to  the 
value  of 


pgv  + / (pUK21  + pWK23)h2dy  . 
o 

determined  retrospectively  from  the  boundary-layer  calculation  at  each 
x-z  station,  whereupon  Equation  (2.74)  reduces  to  a partial-differen- 
tial equation  for  6*.  It  is  integrated  over  the  two-dimensional  domain: 
x,z,  subject  to  the  same  side  boundary  conditions  as  were  imposed  on  the 
principal  governing  equations. 

The  components  of  the  wall  shear  stress,  or  skin  friction,  are 
determined  retrospectively  from  the  boundary-layer  velocity  profiles. 
Specifically: 

T = p (3U/3y)  (2.75) 

w w 

X 

i = p (3W/3y)  , (2.76) 

w w 

2 

where  u is  the  molecular  viscosity  at  the  wall, 
w 
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Separation  of  the  boundary  layer  is  deemed  to  occur  when  the  component 
of  the  wall  shear  stress,  normal  to  the  local  sweep  line,  falls  to  zero. 
This  interpretation  of  separation  is  consistent  with  the  customary 
interpretation  of  oil-flow  visualization  pictures.  Its  usefulness,  as 
a criterion  of  separation,  clearly  deteriorates  as  it  becomes  more  diffi- 
cult to  define  the  sweep  lines  — as  on  a wing  of  high  taper,  for  exam- 
ple. Furthermore,  the  more  three-dimensional  the  flow  becomes,  the  less 
certainty  there  is  that  surface-flow  diagnostics  correlate  with  the 
behavior  of  the  flow  in  the  outer  part  of  the  boundary  layer;  it  is  the 
angle  between  the  externul  streamlines  and  the  body  surface  which  con- 
trols the  departure  of  the  pressure  distribution  and  loading  from  their 
lnvlscld  counterparts.  The  proper  definition  of  separation,  from  an 
engineering  standpoint,  should  therefore  focus  on  abrupt  changes  of  this 
streamline  angle,  rather  than  on  conditions  near  the  wall.  Unfortunately 
there  is,  as  yet,  no  agreed-upon  criterion  of  separation  which  reflects 
this  more  realistic  focus;  the  conventional  wall-shear-stress  criterion 
has  therefore  been  adopted  as  a surrogate. 
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SECTION  III 

GEOMETRY  REPRESENTATION 

Several  approaches  are  available  for  representing  the  body 
surface  In  terms  of  general  curvilinear  coordinate  systems.  The  pres- 
ent authors  take  the  view,  however,  that  orthogonal  systems  are  best 
silted  to  boundary-layer  calculations  because  the  governing  equations, 
so  formulated,  preserve  the  physical  clarity  which  Is  essential  to 
proper  description  and  modeling  of  the  flow.  The  system  chosen  here  is 
based  on  local  development  of  the  surface  as  defined  by  a set  of  air- 
frame reference  coordinates  X,  Y,  Z.  A conical  development  is  performed 
at  each  surface  mesh  point  as  depicted  in  Figure  5.  The  figure  shows 
the  arrangement  of  surface  mesh  points  on  an  arbitrary  body  in  the  neigh- 
borhood of  point  A,  about  which  a local  development  is  to  be  made.  The 
developments  are  performed  row-wise,  marching  in  the  general  direction 
of  the  freestream.  This  procedure  will  be  discussed  later.  Of  interest 
for  the  development  at  A are  the  points  B through  F,  of  which  points  C 
through  E lie  generally  "upstream"  of  A. 

The  lengths  AB,  AG,  Al),  AE,  and  AF  are  preserved  in  the  devel- 
opment. In  addition  the  spatial  angles  between  adjacent  pain-  of 
directed  line  segments  are  preserved,  so  that  local  surface  areas  are 
approximately  preserved  as  well.  Figure  6(a)  depicts  the  development  of 
the  region  shown  in  Figure  5.  The  conical  development  generates  radii 
of  curvature  ^ and  r,,  corresponding  to  the  regions  on  the  left  and  right 
hand  sides  of  point  A.  The  appropriate  coordinate  system  is  then  polar, 
with  the  neighboring  mesh  points  located  via  a radius,  r,  and  an  angle,  0. 
The  origin  of  coordinates  is  fixed  by  a line  0 m constant  passing  through 
the  points  B and  C for  the  left  side  and  F.  and  F for  the  right  side.  The 
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two  coordinate  systems  are  thus  arranged  so  that  point  A lies  at  (r^,  0) 
and  (r^,  0)  respectively.  This  arrangement  has  the  desirable  feature 
that  vector  quantities  at  point  A do  not  require  rotation  in  passing 
from  one  coordinate  system  to  the  other.  Local  vectors  AB,  AC,  and  so 
on,  are  generated  from  the  input  ordinate  data.  The  appropriate  dis- 
tances | AB | through  |AF|  are  generated  by  inner  products  such  as 

| AB | - /AB  • AB  , (3.1) 


The  local  coordinates  (r^O.^)  and  (r2»  e2)  are  then  8®nerated  by  extend- 
ing lines  BC  and  EF  to  intersection  with  line  DA  extended.  Figure  5(b) 
depicts  a local  development  in  the  nomenclature  appropriate  to  the  fin- 
ite differencing  scheme  used  in  the  governing  equations.  The  quantity 
\ on  the  left  and  on  the  right,  is  computed  for  each  local  development 
once  the  coordinates  of  points  B and  F are  determined.  The  region  to 
the  left  of  point  A is  taken  to  be  in  the  direction  of  increasing 
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ordinate  z.  All  geometrical  quantities  can  then  be  referred  to  as  "+" 
subscripted  on  the  side  of  increasing  2 and  on  the  decreasing  side. 
In  the  finite  difference  equations  the  various  geometric  coefficients 
may  now  be  identified  from  the  local  development: 


hx  « | DA | 

(3.6) 

h + - r 0X 

B 

(3.7) 

h“  = r2e2 

(3.8) 

F 

K13  “ 1/rl 

(3.9) 

K13  * 1/r2  ’ 

(3.10) 

As  the  implicit  finite  difference  equations  are  developed  along 

the  span  at  each  downstream  station,  it  becomes  necessary  to  form  first 

z-derivatives  of  the  velocity  components  at  every  point  along  the  span. 

The  situation  is  depicted  in  terms  of  the  local  z-differencing  index  n, 

in  Figure  7.  The  figure  depicts  the  development  at  point  n of  the  + and 

- regions.  Also  shown  are  the  appropriate  developments  for  points  (n-1) 

and  (n+1).  It  is  clear  from  the  figure  that,  for  example,  the  components 

U .,  W ,,  computed  in  the  frame  of  reference  at  station  (n-1)  can  not 
n-1  n-1 

simply  be  transferred  to  the  frame  of  reference  at  position  n without 
taking  account  of  the  relative  rotation  between  the  two  frames  of  refer- 
ence. Thus  in  forming  z-dif f erences  of  velocity  at  station  n involving 
velocity  at  station  (n-1)  the  angle. 


r._  :: ___  '-JI 
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(3.11) 


n-1 


2n 


must  be  used  to  rotate  the  velocity  vector  from  station  (n-1).  Scalar 

5 T 

quantities  are  not  affected  so  that  the  four-vector  F , - (U.  W.  T.  q ) 

n-1  ’ * * M n-1 

occurring  in  the  difference  equations  is  given  in  the  local  frame  by 


F . = R • F , 
n-1  OT  n-1 

where 


CCS 

Vi 

- sin  y . 
n-1 

0 

0 

sin 

'n-1 

cos  Y , 

n-1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

(3.12) 


(3.13) 


A similar  rotation  matrix  applies  when  the  z-derivative  at  n involves 
values  at  (n+1) . 

A similar  situation  exits  when  transferring  vectors  in  the 
x-direction.  Figure  8a  shows  a different  view  of  the  development  region 
around  point  A.  Recognizing  that  the  point  D was,  for  the  previous  x- 
step,  the  point  about  which  development  was  performed,  it  is  seen  that 
a streamwlse  rotation  may  exist  between  the  two  coordinate  systems. 

Since  the  only  quantity  needed  to  transfer  from  point  D to  point  A is 
the  relative  rotation  of  the  two  systems,  it  is  sufficient  to  construct 
a tangent  plant  at  point  D and  to  perform  rotations  into  this  plane. 

To  fix  matters,  the  rotation  of  the  planes  containing  DD'  and  DA  is 
performed  about  the  local  sweep  line  in  the  tangent  plane.  The  result  is 
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depicted  in  Figure  8b.  The  heavy  lines  depict  the  vector  distances  as 

seen  normal  to  the  tangent  plane.  The  dashed  lines  are  the  result  of 

rotation  into  the  tangent  plane.  The  angles  y,  and  y are  then  obtained 

a a 

by  forming  a vector  along  the  sweep  line,  CE,  and  computing  the  dot 


products. 


Fd  • CE  = | D’  D | l CE  | cos(|  + y ) 


AI)  • CE  = | AD  | | CE  | cos  - y ) 


(3.14) 


(3.15) 


The  required  rotation  angle  is  then 


Yb  “ Ya  ' Yd  * 


when  the  streamwisc  slope  of  the  tangent  plane  is  positive,  and 


(3.16) 


yb  ‘ " Ya  * 


(3.17) 


when  the  streamwise  slope  of  the  tangent  plane  is  negative. 

Then  in  order  to  transfer  F at  l-l  to  the  coordinate  system  at  t, 
the  following  rotation  is  performed. 


F » r • p 
Vl  OB  rt-l 


(3.18) 


where 


cos  y.  sin  y.  0 
b b 


sin  y.  c°s  Y.  0 
b b 


0 1 


0 0 


(3.19) 


All  the  geometric  quantities  discussed  in  this  section  are  inde- 
pendent of  the  solution  procedure  once  the  computing  stations  have  been 


settled  upon.  Thus  for  a particular  wing  the  geometry  need  only  be 
computed  once,  and  then  stored  appropriately.  The  data  are  called  back 
one  spanwise  row  at  a time  as  the  solution  progresses  over  the  chord. 
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SECTION  IV 
NUMERICAL  SCHEME 


> 


The  governing  equations  are  solved,  by  an  implicit  finite-dif- 
ference method,  erected  upon  the  surface  on  which  the  boundary  layer  is 
developing.  A forward-marching  procedure  is  adopted,  advancing  in  some 
direction:  x,  in  which  the  mean  velocity  component  is  everywhere  non- 
negative. 

1.  ITERATION  PROCEDURE 

The  matrix-vector  equation  (Equation  2.32),  representing  the  four 
principal  governing  equations,  is  nonlinear  because  the  coefficient 
matrices  depend  on  F.  Accordingly,  the  customary  local  linearization 
approach  is  employed,  involving  an  iteration  procedure  which  successively 
updates  coefficient  matrices  until  they  are  consistent  with  the  converged 
values  of  F at  the  station  where  the  solution  is  being  sought.  Flow 
charts  of  the  iteration  procedure  are  shown  in  the  Appendix. 

The  external  flow  is  solved  first,  to  provide  converged  values 
of  Ue,  W^,  in  terms  of  the  local  pressure  coefficient:  C^.  For  this 
purpose  Equation  (2.32),  which  degenerates  to  a partial-differential 
equation  in  only  two  independent  variables:  x,z,  is  integrated  implicitly 
along  a line  of  roughly  constant  x.  The  values  of  U^,  produced  are 
then  used  to  3cale  the  first-guess  solution  within  the  boundary  layer. 

Within  the  boundary  layer.  Equation  (2.32)  is  integrated  by  means 
of  an  alternating-direction  implicit  (ADI)  method,  covering  each  suc- 
cessive cross-plane  of  roughly  constant  x.  Hie  ADI  method  is  embedded 
within  an  iteration  loop,  again,  for  the  purpose  of  updating  the  coeffi- 
cient matrices  to  a status  consistent  with  the  final  solution:  F.  The 
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continuity  equation  is  also  integrated,  within  this  iteration  loop,  to 
provide  current  values  of  the  effective  normal  mass  flux:  pV  + p'v. 

2.  DIFFERENCE  APPROXIMATIONS 

The  integration  mesh  is  configured  such  that  the  forward-march- 
ing procedure  advances  along  lines  of  constant  y and  z.  Accordingly, 
if  i,  m,  n denote  indexing  in  the  x-,  y-,  z-directions , respectively 
(Figure  9),  x-derivati ves  can  be  expressed  simply  as 


<!£><*> 


J-  {F(e)  _ F(e-1)) 

Ax  m,n  m,n 


(4.1) 


where  Ax  is  the  increment  in  x: 


Ax 


mtn 


U-D 

m,n 


(4.2) 


The  cross-plane:  <,  does  not  in  general  correspond  to  a plane 
of  constant  x,  and  so  the  first-order  backward-difference  approxima- 
tion to  3F/3z  has  to  be  written 


(if)(°  „ -L.  { (l-x-)F(i)  + X'F(i”1)  - FU)  ) , (4.3) 

3*  m,n  Az“  m*n  m»n  m,n-l 


where 


and 


X~  - — (x^  - x<*>  } 

Ax  ' m,n  m.n-l1  ’ 


Ax'  - z<‘>  - z<‘>  . 
m,n  m,n-l 


(4.4) 

(4.5) 


Correspondingly,  the  first-order  forward-difference  approximation  to 
3F/3z  is  written 
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<t§>(,)  . -L.  „<*)  + (l . ,V‘>  . , 

32  m.n  *.+  ®,n+l  m,n  m,n 


“•»  Az 


where 


A+  - <x(4)  - xU)  } , 

Ax  m,n  m,n+l 


(D  AD 

m,n+l  " ra.n  * 


(4.6) 


(4.7) 


(4.8) 


In  Equation  (2.32),  3F/3z  is  approximated  by  a backward,  or  forward, 
difference  according  to  the  sign  of  A^:  backward  if  A^  > 0 and  forward 
if  Aj  < 0.  Thus  the  derivative  is  always  formed  in  the  upwind  direc- 
tion, so  as  to  preserve  convective  stability.  Elsewhere,  the  arithmetic 
mean  of  the  backward  and  forward  difference  is  taken;  it  will  be  noted 
that  this  procedure  is  not  necessarily  equivalent  to  a central  differ- 
ence because  Az  is  not  necessarily  equal  to  Az+. 

In  order  to  express  the  first  derivative  with  respect  to  y,  in 
difference  form,  it  is  convenient  to  write  the  term  A^3F/3y,  in  Equation 
(2.32)  as 


A „ fA<l>  + a(2),  3F 

A4  3y  ”4  A4  1 3y  ’ 


(4.9) 


where  A^^  contains  only  diagonal  elements  and  A^^  only  off-diagonal 


elements.  Then: 


. ,3F ,<*>  A4  ...  . t,AD  ,.pU) 

A.(t~)  * T 1(4’  + DF^,  - 2<J>F 

4 3y  m,n  2 m+l,n  m,n 


+ <*-l>^l.nl  +T  'F^,n-Fm-i,nl- 


(4.10) 


The  quantity  ♦ , in  Equation  (4.10)  is  a diagonal  matrix,  equivalent  to 
the  parameter  a employed  in  Reference  [16],  which  controls  the  level  of 
artificial  viscosity  to  be  introduced  via  the  difference  equations.  The 


elements  of  <^,  lie  between  -1  and  1,  and  have  the  same  sign  as  the 
corresponding  terms  of  The  difference  equations  become  uncondi- 
tionally stable  if  = 1,  but  are  only  conditionally  stable  if 

< 1.  It  will  be  noted  that  no  artificial  viscosity  is  introduced 
via  the  off-diagonal  terms  of  ; this  is  appropriate  since  any  artifi- 
cial viscosity  so  introduced  would  relate  one  element  of  3F/3x  to 
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another  element  of  3 F/3y  , with  unpredictable  results. 

Second  derivatives  with  respect  to  y are  expressed  in  the  usual 


~ ' — " 2p(i)  + f(*i  (A. 11) 

.2  .2  ro+l,n  m,n  m-l,n 

3y  m,n  (Ay) 


where 


(O  (t)  (l)  (l) 

,y  ym+l,n  ym,n  ym,n  ym-l,n* 


(4.12) 


Upon  introduction  of  the  various  difference  approximations  into 
Equation  (2.32),  the  following  difference  equation  is  obtained,  relating 
to  the  node  point  f,  m,  n (see  Figure  9): 


_(*)  n F(0  . d -(O 

BlFnH-l,n  Tm,n+1  3 m,n 


+ B,F(*j  + BF(l)  . - B F(l_1),  (4.13) 

4 m-l,n  5 m,n-l  6 m,n 


where  the  B's  are  square  matrices  representing  linear  combinations  of 
the  A's  in  Equation  (2.32). 

3.  ALTERNATING-DIRECTION  SCHEME 

The  forward-marching  procedure  advances  in  the  positive  x-direc- 
tion,  from  a cross-plane:  l-l , where  the  solution  is  assumed  to  be 
known,  to  a cross-plane:  £,  at  which  a new  solution  is  sought.  The 
cross-plane  is  scanned  alternately  in  the  n-,  and  m-directions  (Figures 
10,  11),  converting  Equation  (4.13),  respectively,  into  the  successive 
forms : 


1 F^1 

+ B_F^  + 

B F^| 

1 m+l,n 

3 m,n 

4 m-l,n 

FU) 

+ B_F^  + 

B,FU) 

2 m,n+l 

3 m,n 

5 m,n-l 

(4.14) 

(4.15) 


in  which  C^,  C^  contain  the  passive  terms  originating  on  the  left-hand 
side  of  Equation  (4.13).  During  the  scan  in  the  n-direction.  Equation 
(4.14)  represents  a recursion  relationship  for  1 £ m M,  and,  during 
the  scan  in  the  m-direction,  Equation  (4.15)  represents  one  for  1 ± n N. 
In  either  case  the  sequence  of  Equations  (4.14,  15)  may  be  written  in 
the  form 


PQ  - R , (4.16) 

where  Q is  a vector  of  order  M or  N,  whose  elements  are  the  vectors  F, 
R is  a vector  of  similar  order  containing  the  elements  C^  or  C2,  and  P 
is  a tridiagonal  square  matrix,  of  order  M or  N,  containing  the  B's 
from  Equation  (4.14)  or  (4.15). 


- 1 

Che  solution  of  Equation  (4.16)  is  formally  determined  as 


Q - P_1  R . (4.17) 


L. 

The  solution  is  readily  calculated  using  the  extended  Choleskl 
method  (Reference  [17])  which  involves  forward  elimination  of  the  lower 


diagonal  elements  of  P,  and  backward  substitution  to  determine  the  ele- 


ments of  Q.  The  coefficient  matrices  P,  R are  updated  in  successive 
iterations  and  the  solution  process  continued  until  convergence  is 
obtained. 
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SECTION  V 

ADDITIONAL  DETAILS  OF  THE  METHOD 


1.  EXTERNAL  FLOW 

Ic  Is  assumed  that  the  external  flow  conditions  are  specified 

in  terms  of  the  local  pressure  coefficient,  C , together  with  condi- 

P 

tions  in  the  undisturbed  free  stream.  The  local  static  pressure, 
static  temperature,  Mach  number,  and  resultant  velocity  in  the  exter- 
nal flow  can  be  determined,  directly,  as  follows: 


-E_  = 
Pcc 


1 + 


T 

e 

T 


III 
(JL)  Y 


e 


2 

Y-l 


{(1  + 


(5.1) 


(5.2) 


(5.3) 


Q M T 1/2 
e e e 


(5.4) 


The  velocity  components,  Ue>  W^,  cannot  be  determined  directly,  but  must 
be  found  by  solving  the  Euler  equations  (the  momentum  equations  without 
the  viscous  and  turbulent  shear-stress  terms)  subject  to  the  appropriate 
boundary  conditions.  These  equations  are  solved  implicitly,  in  the 
two-dimensional  domain:  x,  z,  using  the  same  difference  approximations 
as  are  employed  for  the  full  boundary- layer  equations. 

Truncation  errors  associated  with  these  difference  approxima- 
tions will  of  course,  arise  and,  as  a result,  the  values  of  Ug  and  Wg 
derived  will  not  generally  be  consistent  with  the  value  of  Qg  determined 
from  Equation  (5.4).  The  truncation  errors  can  however,  be  reduced  for 
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certain  types  of  flow,  by  suitable  handling  of  the  pressure  gradients: 
3p/3x,  3p/3z.  In  particular,  for  an  infinite  yawed  wing,  a transforma- 
tion of  the  form: 


C2  - cos2A - 

PT 


(y-i)m 


yM2  ^ 7 

(1  + tV  “ 1 (*  (5‘5> 


will  enable  the  exact  values  of  11  , W , to  be  obtained  using  a two- 

e e 

point  backward  difference  for  3C  /3x: 

P 

3x  Ax  ) 2 p ) pT  n PT  n ( 


and  a three-point  difference  for  3C  /3z: 
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3C  yM2  \ i f|v  . ,+ 
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The  transformed  pressure  coefficient:  C , will  be  recognized  as  being 

PT 

Identical  to  the  velocity  component  normal  to  the  sweep  lines,  and  the 
quantity 


1/y 

(1  ♦ — cp) 


in  Equations  (5.6,7)  as  being  Identical  to  the  local  density  ratlo:p  /p^. 
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In  Che  present  formulation,  the  above  transformation  is  used 
for  the  finite  swept  wing  also.  In  this  more  general  case  truncation 
errors  will  arise  in  the  solution  of  the  Euler  equations  but,  in  most 
cases,  they  will  be  smaller  than  the  errors  that  would  arise  if 
3C^/3x,  3Cp/3z  were  determined  by  applying  the  difference  formulae  to 
the  untransformed  pressure  coefficients.  It  should  be  stressed  however, 
that  no  approxlmat Ion  is  Introduced  into  the  Euler  equations  other  than 
the  difference  approximations  to  the  derivatives. 

2.  DISTRIBUTION  OF  NODE  POINTS  THROUGH  THE  BOUNDARY  LAYER 

The  thickness  of  the  integration  domain  is  set  equal  to  ftfl , where 
6 is  the  boundary  layer  thickness  and  i)  is  a constant  in  the  neighbor- 
hood of  1.23.  Thus,  if  the  physical  distance  from  the  wall  is  denoted 
by  y,  then 

0 < y < ftS  as  0 < m < M , 

where  M is  the  number  of  node  points  distributed  through  the  domain. 

The  transformed  space  variable,  y,  is  defined  by 

y « rnAy  , (3.8) 

where  Ay  is  any  convenient  constant,  ard  hence  the  relationship  between 
y and  y describes  the  manner  in  which  the  node  points  are  distributed. 

For  laminar  flow,  a uniform  distribution  is  appropriate,  where- 
upon 

y ■ ft6m/M  , (5.9) 

h.,  “ 3y/3m  * ftd/M  , (5.10) 

A. 

K - 0.  (5.11) 
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For  turbulent  flow,  a non-uniform  distribution  is  more  appropriate 
which  provides  Increased  point-density  in  the  region  of  large  gra- 
dients close  to  the  wall.  A distribution  that  respects  the  various 


regions  of  a turbulent  boundary  layer:  viscous  sublayer,  logarithmic 
region,  and  "wake  region"  is 

y “ -(r£V)M  U^m-M)  + AM)  , (5.12) 

in  which  <p  m 4>(m)  is  such  that  $ becomes  large  as  m+M.  The  distribu- 
tion thus  reduces  to  the  uniform  one 

y - ~ (c1(m-M)  + AM}  , (5.13) 

in  the  outer  part  of  the  boundary  layer.  In  the  inner  part,  the  dis- 
tribution is  controlled  by  the  function  $ which  takes  the  form 

cm 

$ - c2m  + c^(e  4 - 1)  . (5.14) 

Here,  c^  is  taken  to  be  proportional  to  the  additive  constant  in  the 
Law  of  the  Wall,  and  c2  determines  the  distance  from  the  wall  at  which 
it  is  desired  to  place  the  first  node  point;  c^  is  selected  such  that 
y is  non-singular  in  the  interval  0 <_  m M. 


The  values  of  h2  and  K22,  corresponding  to  Equation  (5.12)  can 
be  determined  by  differentiation  with  respect  to  m: 
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SECTION  VI 

EVALUATION  OF  THE  METHOD 


Time  and  funding  constraints  have  placed  severe  limitations  on 
the  scope  of  the  validation  procedures  which  could  be  performed  on  the 
method.  The  following  represent  only  a cursory  set  of  checks  to  deter- 
mine whether  gross  errors  were  likely  to  arise  when  the  method  was 
used.  A more  comprehensive  process  of  validation  is  clearly  needed,  to 
search  for  more  subtle  errors  or  inaccuracies,  involving  checks  of 
internal  consistency  as  well  as  comparisons  with  experiment  and  exact 
theory.  These  comparisons  would  also  form  a basis  for  "fine  tuning"  of 
the  method  (adjustment  of  the  empirical  functions,  etc.)  which,  so  far 
has  not  been  done. 

A number  of  numerical  convergence  checks  were  made  using,  as  a 
simple  test  case,  two-dimensional,  incompressible  laminar  flow  over  a 
flat  plate.  Figure  12  shows  the  effect  of  varying  the  number  of  node 
points  distributed  through  the  boundary  layer.  It  is  evident  that  20 
points  is  sufficient  to  adequately  recover  the  exact  Blasius  solution 
and  that  further  refinement  has  no  significant  effect.  The  results  shown 
were  obtained  using  7 streamwise  stations  (Including  the  leading  edge 
which  was  represented  as  a stagnation  point  from  which  the  fluid  accel- 
erated linearly  to  free-stream  velocity  in  a distance  equal  to  0.005  times 
the  length  of  the  plate) . A coarser  description  involving  only  5 stream- 
wise  stations  resulted  in  a 4%  error  in  wall  velocity  gradient  compared 
to  an  error  of  less  than  1%  with  7 stations.  A finer  description, 
involving  14  streamwise  stations,  reduced  this  error  to  less  than  0.05Z. 

Figure  13a  shows  the  convergence  of  the  iteration  process  for 
the  case  involving  7 streamwise  stations  and  20  node  points  across  the 


A. 
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boundary  layer.  At  this  particular  x-station  (the  7th)  the  calculation 

-4 

converged  to  a tolerance  of  10  in  6 Iterations.  Here,  the  tolerance 

is  defined  as  the  maximum  permissible  variation  between  successive 

iterations,  expressed  as  a fraction  of  Ug.  (Corresponding  tolerances 

~~2 

are  placed,  where  applicable,  on  W,  T,  and  q ).  Relaxation  of  the  tol- 
erance results  in  progressively  larger  errors  in  wall  velocity  gradient 
(Figure  13b);  with  a tolerance  of  0.003  this  error  is  about  8X. 

Figure  14a  shows  the  predicted  effect  of  compressibility  on  the 
two-dimensional  laminar  boundary  layer.  The  calculations  indicate  a 
reduction  of  about  162  in  normalized  wall  shear  stress  at  a Mach  number 
of  3.  This  order  of  reduction  is  close  to  the  figure  given  by  Young's 
formula  (18].  Figure  14b  shows  the  computed  temperature  profile  through 
the  laminar  boundary  layer  at  Mach  numbers  of  2 and  3.  The  profiles  show 
the  enthalpy  is  being  properly  conserved  in  the  outer  part  of  the  bound- 
ary layer.  The  loss  of  enthalpy  near  the  wall  is  of  course  to  be 
expected  because  the  recovery  factor  is  less  than  one. 

Figure  15  shows  some  calculated  results  for  incompressible  lami- 
nar flow  over  an  infinite  45°-swept  wedge.  The  external  velocity  compo- 
nent: U , normal  to  the  leading  edge,  is  assumed  to  increase  in  propor- 
e 

n 

tion  to  the  normal  distance:  x , while  the  external  component:  W , 

n e 

P 

parallel  to  the  leading  edge  is  constant.  The  computed  velocity  pro- 
files are  expressed  in  similarity  form  and  are  compared  with  the  exact 
solution  of  Howarth  and  Cooke  [19,20].  It  should  be  noted  that  the  cal- 
culation proceeded,  not  along  the  normal  to  the  leading  edge,  but  in  the 
direction  of  the  undisturbed  stream.  The  velocity  profiles:  U(y),  W(y), 
appearing  in  the  calculation  did  not,  therefore,  reach  an  asymptotic  form 
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but  varied  with  x.  The  profiles  plotted  in  Figure  IS  were  obtained  by 
resolution. 


Some  results  for  incompressible  turbulent  flow  are  presented 
in  Figures  16  through  19.  Figure  16  shows  an  inner-law  plot  of  the 
velocity  profiles  in  a constant-pressure  turbulent  boundary  layer  at  4 
Reynolds  numbers.  It  will  be  noted  that  the  viscous  sublayer  and  blend- 
ing region  are  well  represented,  and  that  the  familiar  logarithmic  region 
is  obtained  at  high  Reynolds  numbers.  The  slope  of  the  logarithmic  por- 
tion is  consistent  with  the  value  of  von  Karman's  constant  used  here 
(0.418  );  the  additive  constant  emerges  as  2.5,  compared  with  Patel's  [4] 
value  of  2.3. 

The  predicted  wall  shear  stress  is  expressed  as  a function  of 
Reynolds  number  in  Figure  17.  The  trend  is  consistent  with  that  given  by 
Rotta  [21],  Coles  [22],  and  Smith  and  Walker  [23],  although  the  predicted 
values  themselves  appear  to  be  slightly  low. 

Figure  18  shows  the  computed  mean  velocity  profile  at  a Reynolds 

4 

number  (Rfi*)  of  10  . Two  sets  of  results  are  plotted,  corresponding  to  20 

and  30  node  points  across  the  boundary  layer,  and  there  is  no  significant 

difference  between  them.  The  calculations  are  in  good  agreement  with  the 

experimental  data  of  Klebanoff  [24]  and  Smith  and  Walker  [23].  The  corre- 

~ 2 

sponding  profiles  of  turbulence  intensity:  q , are  shown  in  Figure  19. 

It  will  be  noted  that  q Increases  from  zero,  at  the  wall,  to  a maximum 

(for  constant  - pressure  flow)  which  is  located  near  the  edge  of  the  vis- 

2 2 

cous  sublayer.  The  maximum  value:  q /Ug  “ 0.0085,  is  consistent  with 

— ~2 

the  predicted  wall  shear  stress  of  0.0013,  assuming  a value  of  aj  («-uv/q  ) 
equal  to  0.15.  The  calculated  results  are  compared  with  the  data  of 
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Klebanoff  [24].  Experimental  values  of  qZ  are  shown  and  also  the  meas- 
ured values  of  -uv  divided,  here,  by  a^  The  two  sets  of  experimental 
data  are  in  close  agreement  with  each  other  over  most  of  the  boundary 

layer,  but  diverge  close  to  the  wall  because  of  the  "inactive"  component 

~2 

of  the  turbulence  included  in  q . The  calculated  results  are  in  closer 

agreement  with  the  values  of  -uv/a^  in  this  region,  confirming  that  the 
~2 

quantity  q considered  in  the  theory  is  to  be  identified  with  the  active 
component  of  the  turbulence  which  is  responsible  for  momentum  transport. 
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SECTION  VII 

DEMONSTRATION  TEST  CASE 

The  test  case  selected  for  demonstration  purposes  corresponds 
to  the  upper  surface  of  the  right  wing  of  the  NASA/Air  Force  TACT 
(modified  F 111)  airplane.  The  wing  geometry  and  flight  pressure  dis- 
tributions have  been  supplied  by  NASA  Flight  Research  Center,  Edwards, 
California.  The  particular  case,  for  which  a boundary  layer  calcula- 
tion has  been  performed  corresponds  to  a nominal  leading  edge  sweep  angle 

of  26°,  and  an  airplane  angle  of  attack  of  9.01°.  The  Mach  number  is 

2 

0.84,  and  the  static  pressure  is  903.08  lb/ft  . 

The  wing  geometry  is  shown  in  Figure  20.  The  pressures  were 
measured,  in  flight,  along  four  spanwise  stations  designated  Rows  A 
through  D,  in  which  Row  A is  near  the  tip,  and  Row  D near  the  junction 
with  the  inner-wing  glove.  The  instrumented  portion  of  the  wing  was 
clean,  and  was  uniformly  tapered.  The  pressure  distributions,  along  the 
four  spanwise  stations  are  shown  in  Figure  21,  and  it  will  be  noted  that 
a shock  is  present  over  most  of  the  span.  The  shock  appears  to  be  simple 
at  Rows  B and  C,  but  becomes  bifurcated  at  Row  A,  near  the  wing  tip. 

The  calculation  was  performed  with  15  streamwise  stations  and  four 
spanwise  stations  (corresponding  to  Rows  D,  B,  C,  A).  The  end  stations 
(D,A)  were  treated  as  being  portions  of  locally-infinite-swept  wings. 

In  the  interior  of  the  integration  domain  (Rows  C,B)  the  flow  was  treated 
as  being  fully  three-dimensional.  Transition  was  assumed  to  occur  at 
11%  of  local  streamwise  chord,  on  Rows  D,  C,  and  at  6%  of  local  chord 
on  Rows  B,A.  These  stations  correspond  to  the  start  of  the  adverse 
streamwise  pressure  gradient. 
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and 


Figures  22  and  23  present  the  distributions  of  6*,  x , and 

x 

t versus  local  chord  at  the  four  spanwise  stations.  Figure  24  presents 
z 

velocity  profiles  as  computed  at  three  representative  points  on  the  plan- 
form.  Separation  occurred  in  a small  region,  depicted  approximately  in 
Figure  24,  and  the  separation  control  process  was  used  to  recover  the 
boundary  layer  calculation.  Evidence  of  approaching  separation  can  be 
seen  in  Figures  22  and  23  through  the  behavior  of  6*  and  t 
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SECTION  VIII 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  method  described  herein  builds  upon  experience  in  calculating 
three-dimensional  boundary  layers  extending  over  some  eight  years.  This 
experience  has  been  reflected  in  the  retention  of  many  of  the  more  firmly- 
based  features  of  the  earlier  methods,  and  these  features  have  been  incor- 
porated, together  with  many  new  ones,  into  an  entirely  new  program  struc- 
ture. The  new  structure  offers  immediate  computational  advantages  with 
respect  to  the  ease  of  setting  up  a calculation  and  the  economy  of  the 
calculation  itself.  It  also  offers  considerable  flexibility  which  will 
permit  extensions  and  refinements  of  the  method  as  the  need  for  them 
arises.  In  that  regard  the  method  may  be  considered  to  be  a starting 
point  for  further  work,  rather  than  a completed  product,  even  though  it 
offers  the  user  a comprehensive  tool  which  should  more  than  adequately 
perform  the  tasks  for  which  it  is  currently  needed.  It  will  clearly  take 
time  to  properly  evaluate  the  new  method,  and  to  exploit  its  full  poten- 
tial. 

Central  to  the  new  program  structure  is  the  alternating-direction 
implicit  (ADI)  numerical  scheme.  We  believe  that  the  adoption  of  this 
type  of  numerical  scheme  represents  a definitive  solution  of  several  of 
the  problems  which  have  plagued  the  development  of  implicit  methods  for 
calculating  three-dimensional  boundary  layers.  We  also  believe  that  it 
represents  the  best  starting  point  for  several  potential  extensions  of 
the  method,  involving,  for  instance,  the  inclusion  of  higher-order  terms 
in  the  equations  of  the  mean  motion  and/or  the  turbulence  model  equations. 

A number  of  further  refinements  to  the  numerical  scheme  would, 
however,  be  worthwhile  even  for  flows  which  satisfy  the  first-order 
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boundary-layer  equations.  In  particular.  It  would  be  advantageous  to  up- 
grade the  schene  to  second-order  accuracy  with  respect  to  the  x-  and  z- 
derlvatlves.  Also,  improvements  could  be  made  to  accelerate  the  conver- 
gence of  the  Iteration  process.  Schemes  which  claim  to  eliminate  the 
need  for  the  iteration  process,  altogether,  should  be  evaluated  too;  but 
the  present  authors  are  skeptical  about  the  desirability  of  such  schemes 
if  they  place  severe  restrictions  on  the  step-size  in  the  direction  of 
march.  If  this  is,  indeed,  the  case,  the  iteration  process  has  simply 
been  traded  for  a much  larger  number  of  x-steps. 

One  of  the  features  of  the  earlier  methods,  which  has  been 
retained,  is  the  one-equation  turbulence  model.  We  have  consistently 
argued  the  superiority  of  such  a model,  over  the  "zero-equation"  eddy- 
viscosity  or  mixing- length  approaches,  on  the  grounds  of  its  greater 
universality.  On  the  other  hand,  this  universality  could  be  improved 
still  further  by  the  Incorporation  of  at  least  one  more  differential  equa- 
tion: for  the  dissipation  rate  or,  equivalently,  the  turbulence  length 
scale.  We  recommend  that  this  refinement  be  made,  possibly  in  an  experi- 
mental method,  to  enable  the  user  to  check  the  validity  of  the  existing 
model  under  specified  flow  conditions  where  its  adequacy  might  be  ques- 
tionable. 

To  the  extent  that  the  turbulence  model  is  closely  similar  to 
that  used  in  the  earlier  methods,  it  would  be  expected  that  the  degree  of 
correlation  between  the  predicted  mean-flow  quantities  and  experiment, 
would  be  comparable  to  that  previously  reported  for  turbulent  flow  [1,3,4]. 
A wider  range  of  measurements  is  now  available,  however,  and  it  would  seem 
highly  desirable  to  make  comparisons  with  these  recently  acquired  data. 


The  comparisons  would  serve  as  a further  check  on  the  reliability  of 
the  method,  and  might  also  indicate  the  areas  where  "fine-tuning"  of 
the  empirical  functions  would  produce  even  better  agreement  with  experi- 
ment. It  is  worth  noting  that  the  assumed  relationship  between  dissipa- 
tion length  and  position  through  the  boundary  layer,  and  the  diffusion 
function  assumed  herein,  are  identical  to  their  counterparts  in  the 
earlier  method  [3,4]. 

Some  improvements  can  still  be  made  to  the  geometry-representa- 
tion scheme.  The  conical-projection  approach  is  probably  adequate  for 
most  purposes,  but  it  could  be  upgraded  to  a full  doubly-curved  repre- 
sentation, with  correspondingly  greater  generality,  without  difficulty. 

The  governing  equations  are  already  cast  in  terms  of  two  geodesic  curva- 

| I 

, tures.  In  this  regard,  however,  we  feel  that  there  is  a need  to  accumu- 

late user  experience  in  the  description  of  typical  airplane  geometries, 
prior  to  the  decision  to  replace  the  existing  scheme. 

One  of  the  major  sources  of  uncertainty,  concerning  a boundary- 
layer  calculation  extending  only  over  the  wings  of  an  airplane  configu- 
ration, has  to  do  with  the  specification  of  realistic  boundary  conditions 
at  the  wing  root.  Fluid  spilling  off  the  fuselage,  onto  the  wing,  convects 
information  which  is  relevant  to  the  boundary- layer  development  but  which 
the  user  is  unable  to  specify.  Ideally,  the  calculation  should  extend 
over  the  whole  of  the  wetted  surface  of  the  airplane,  in  which  case  the 
need  for  specifying  boundary  conditions  along  such  arbitrary  interfaces 
would  be  eliminated.  Unfortunately,  such  interfaces  often  correspond, 
as  does  the  wing-fuselage  junction,  to  regions  where  the  first-order 
boundary-layer  equations  are  invalid,  and  where  a more  sophisticated 

t ■ 
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calculation  would  be  required.  Be  that  as  it  may,  the  present  method 
could  readily  be  extended  to  a wing-fuselage  combination  on  the  assump- 
tion that  the  corner  flow  is  amenable  to  first-order  treatment.  The 
extended  method  would,  at  least,  take  account  of  continuity  of  mass 
across  the  interface,  and  would  partially  represent  the  momentum 
exchange  too.  We  strongly  recommend  that  such  an  extension  of  the  method 
be  funded,  at  an  early  date,  to  provide  an  interim  engineering  tool  for 
airplane  boundary- layer  computations,  and  to  provide  a framework  for 
incorporation  of  the  more  sophisticated  treatment  of  the  corner  flow  when 
the  technology  for  handling  it  becomes  available. 
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forward  view 


igure  1 WING  GEOMETRY,  SHOWING  SEGMENTATION  INTO  PANEL  ELEMENTS 


Figure  3 EMPIRICAL  TURBULENCE  FUNCTIONS 


LATERAL  ROTATIONS  BETWEEN  ADJACENT  ELEMENTS  OF  THE  SAME  ROW 


previous  level 


forward-marching 

direction 


current  level 


Figure  9 FINITE-DIFFERENCE  MOLECULE 
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Figure  11  FINITE-DIFFERENCE  MOLECULE,  SHOWING  ALTERNATING- 

DIRECTION  SCANNING 
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Figure  13  CONVERGENCE  OF  THE  ITERATION  PROCESS 
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DIMENSIONAL,  INCOMPRESSIBLE  TURBULENT  FLOW  OVER  A FLAT  PLATE 
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Figure  20  DEMONSTRATION  TEST  CASE:  WING  GEOMETRY 
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Figure  24  DEMONSTRATION  TEST  CASE:  VELOCITY  PROFILES 
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5.  Description  of  Subroutines 

6.  Program  Flow  Diagrams 


NASH/SCRUGGS  THREE  DIMENSIONAL  BOUNDARY  LAYER  CODE* 


Input  Procedure 

The  input  data  is  contained  in  several  blocks.  All  input  is  read  into 
the  MAIN  program  added  at  Grumman. 


Card  # 

Format 

Field 

Name 

Remarks 

A1 

Literal 

GTITLE1 

1st  card  of  2 card 
title  set 
(STOP  in  the  1st  4 
columns  terminates 
execution) 

A2 

Literal 

- 

GTITLE2 

2nd  card  of  2 card 
title  set 

A3 

2113 

1 

LL 

# of  streamwise 
stations  (=£50) 

2 

MM 

# of  computational 
points  across  bound- 
ary layer  (<21) 

3 

NN 

# of  span  stations 
«30) 

4 

ITN 

Maximum  number  of 
iterations  allowed 
in  numerical  pro- 
cedure at  each  stream- 
wise  row 

A4 

2113 

1 

KBC 

Edge  boundary  conditions 

=1  2-D  flow 

=2  plane-of-symmetry 
flow 

=3  infinite  yawed - 
wing  flow 

2 

IOITER 

Output  of  numerical 
iteration  history 

=0  no  output 

=1  full  output 

3 

IOPRF 

Output  of  velocity 

profiles  at  each 
station 


=0  no  output 
=1  full  output 


* Input  format  revised  slightly  by  Bill  Mason  at  Grumman,  where  this 
section  was  rewritten. 
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BLOCK  B Geometry  Input 


Card  # 

Format 

Field 

Name 

Remarks 

B1 

7F10.0 

1 

xxgeom(n,i) 

X (streamwise)  location 
of  leading  edge  at  span 
station  N 

2 

XXGE0M(N,2) 

Y (vertical)  location 
of  leading  edge  at  span 
station  N 

3 

XXGE0M(N,3) 

Z (span)  location  of 
leading  edge  at  span 
station  N 

4 

xxgeom(n,4) 

Local  chord  at  span 
station  N 

5 

XXGE0M(N,5) 

Wing  incidence  (in 
addition  to  ordinates) 
at  span  station  N 

B 2 

7F10.0 

1-LL 

XYGEOM(Ll,N,l) 

(-)  ordinates  at  span 
c station  N,  read  for 
LI  from  1 to  LL 

NOTE:  Card 

B2  is  read 

seven  values  per 

card  until  LL  values 

are  read. 

B3 

7F10.0 

1-LL 

XYGE0M(L1,N,2) 

(£)  ordinates  at  span 

station  N,  read  for 
LI  from  1 to  LL 

NOTE:  Card  B3  is  read  seven  values  per  card  until  LL  values  are  read. 

NOTE:  BLOCK  B is  repeated  once  for  each  span  station  starting  with  the  most  inboard 
span  station  (NN  times). 


BLOCK  C Flight  Conditions 


Cl 

Literal 

TITLE 

Title  card  describing 
flight  condition  and 
pressure  distribution. 

C2 

7F10.0 

1 

XMINF 

Mach  number 

2 

XLEN 

Reference  length  scale, 
in  feet 

3 

PINF 

Static  pressure,  in  psi 

4 

TINF 

Static  temperature,  in 
degrees  Rankine 

C3 

2113 

1— NN 

LTRN(N) 

Chordwise  location  of 
specified  transition 
at  each  span  station 

NOTE:  Card  C3  is  read  a second  time  if  NN  > 21 


Card  # 


Format 


7F10.0 


Field 


BLOCK  D Pressure  Distribution 


7F10.0 


IPNCH 


KSMTHS 


CNCRIT 


CRTLAM 


CRTTBL 


CP(L,N) 


Remarks 

= 0 do  not  write 
solution  to  file 

= 1 write  solution  to 
file 

= 0 do  not  read  solution 
fran  file 

= 1 read  solution  from 
file 

= 0 do  not  punch  output 
for  plotting 

= 1 punch  output  for 
plotting 

Number  of  smoothings  of 
input  pressure  distribu- 
tion (if  desired) 

Convergence  criterion 
for  numerical  iteration 
(.0005  nominal  value) 

Separation  test  for 
laminar  flow  (.05 
nominal  value) 

Separation  test  for 
turbulent  flow  (.005 
nominal  value) 


Pressure  distribution 
at  constant  span  station, 
L.E.  to  T.E. 


NOTE:  Card  D1  is  read  repeatedly  until  LL  values  are  input,  seven  values 
per  card. 

NOTE:  BLOCK  D is  read  once  for  each  span  station  (NN  times). 

Program  Termination:  The  program  will  automatically  return  to  the  beginning  and 

read  card  Al.  If  the  1st  4 columns  in  Card  A1  contain  the 
clue  work  STOP,  execution  terminates. 
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IBIS  FAGS  IS  BEST  QUALITY  mCTXfiABfeS 

, PRQU  COi  I FURNISHED  TO  DDC  - 

IBM/ CUC  Conversion  v 1 " 

Only  minor  modifications  are  required  to  convert  the  program  from  CDC  to 
1 1 , i-  versa.  The  following  list  provides  the  entire  details  of  tlu 

conversion  procedure. 


I.  Comment  out  the  program  cards  for  IBM. 

II.  Entry  points  must  be  changed.  These  occur  in  the  listed  subroutines: 


SUBROUTINE 


CDC 


IBM 


1.  BOUND 

2.  PRINT 

3.  PRINT 

4.  DISP 

5.  GEOMT 


ENTRY  GUESS 
ENTRY  PRINT.? 
ENTRY  PRINT3 
ENTRY  DISP2 
ENTRY  GE0M2 


ENTRY  GUESS  (L) 
ENTRY  PRINT2(L) 
ENTRY  PRINT3(L) 
ENTRY  DISP2(L) 
ENTRY  GE0M2 ( L ) 


III.  For  those  IBM  installations  that  do  not  have  the  H-Extended  Compiler,  two 


function  subprograms  are  required  in  order  to  supply  ASIN  and  ACOS.  This  is 
done  by  using  the  function  subprograms  to  equate  ASIN  to  ARSIN  and  ACOS  to  ARCOS. 


IBM  JCL 


//G*U»<A001  JOB  (727a»I36.OO50,0M30,O0P0,,2,,62),  ' * AC  KE  NZ  J F •!>75«B722  ' , X 

//  MSGLtVEl«l>Cl.A3S«C*BfcGlON»*?,SfcK,TlMfc*15 

/•MESSAGE  CPU  T I ME  IS  *IM,  aALL  CLOCK  50  mjm,  region  256* 


//NASH  EXEC  E0RTHCLG,PAt?MtEn»T«»llPT«2'  , Rfc'G ION ,FORT«256K , X 

//  BARMiLKfn*'SI/H(H0*,2»5).  HAP,  IFT,  LIST',  x 

//  T T MF , GO* 1 5 

//A  OPT  , S YSL  I N Of'  PSN*A&L0AD8ET,0ISP*(,PASS),  X 

//  UNIT»(SV3DA#S§P*8YSPHINT),.1PACE*(CYI,  («,!)),  X 

//  flCH*(PECEM«»-H,LRECL«H0,BLKSI2fc*3?00) 

//PORT*  SYS  IN  [in  PI3P»SHR,PISN«D77027B, MASON, GACNASH.CAHOS,  X 

//  UNJT«3330, vrt«SER«OSPE02 

//LKfn.SYSIN  00  * 

ENTRY  MAIN 

/A 

//Gn,A TO7P00I  00  dummy 


//GO.ETIUEOOI  DO  UN  I T«S Y$0 A , SPACE ■ (C Yl , ( ? , 1) ) , 0 1 SP« ( # Dfc Lfc TF  ) , X 

//  DC8»(RECAM«v§s,LRf CL«32#BLK3I/F»70a) 

//an, ft»5m»pi  o«  mummy 

//QP.PTjfcPOOl  00  UNlT«8YSOA,SPACr«(TR«, (3, 5)),0ISR»(, DELETE),  X 

//  0CR*(REC»“«^3,LRErL»S?,»L'»3T2E»«6) 

*YO  • 

rC*»r*  •<  tt  *•*•!) 

/• 


1.  List  of  Primary  Coding  Symbc.’Vs 

Indices: 

M index  in  the  y-direction,  MM  upper  limit. 

N index  in  the  z-direction,  NN  upper  limit. 

L index  in  the  x-dlrectlon,  LL  upper  limit. 

IT  iteration  number,  upper  limit  ITN. 

F(M,N, 1)  U(y,z) 

F(M,N, 2)  W(y,z) 

F(M,N, 3)  T(y,z) 

F(M,N,4)  q2(y,z) 

FO(M,N,I)  the  quantities  above,  1-1,4,  for  the  previous  x-station. 

(In  general  an  "0"  appended  will  have  this  meaning.) 

FP(M,N,I)  values  of  F(M,N,I)  from  previous  iteration 

CP(L,N)  C (x,z) 

P 

H1(N)  metric  coefficient  h^. 

H3(M,N, I)  metric  coefficient  (13,  measured  on  the  side  of  increasing 

z (I  - 1)  and  on  the  side  of  decreasing  z (1-2)  at  local 
development  point. 

DISP(N)  displacement  thickness 

LAMINAR(N)  equals  1 for  laminar  flow,  2 for  turbulent  flow. 

LTRN(N)  value  of  L for  which  transition  is  specified,  at  each  N. 

(input) 

KBC  an  integer  fixing  boundary  conditions  in  the  z-direction. 

( input ) 

QE(N)  velocity  at  the  edge  of  the  boundary  layer. 

RHOE(N)  density  at  edge  of  the  boundary  layer. 

R0T1(N,I,J)  rotation  matrix  for  the  side  of  decreasing  z at  the  local 

development  point  (I,J-1,4). 

R0T2(N,I,J)  same  as  R0T1  for  left  hand  side. 

?WP(N)  leading  edge  sweep  angle. 

■Ug(g.I.J)  ha«  toward  rotation  metrta  a the  local  developaient  point 


SWPL(N) 

local  sweep  angle. 

XME(N) 

Mach  number  at  edge  of  the  boundary  layer. 

XK13(N,I) 

the  curvature  coefficient  K..,,  measured  on  the  side  of 
decreasing  z (1*2),  and  on  the  side  of  increasing  z 
(I  * A)  at  the  local  development  point. 

V(M,N) 

pv  + pT*v 

XK2(M,N,I) 

the  curvature  coefficient  K£2  ^or  t*ie  previous 
(1*1)  and  the  current  x-level  (1*2). 

x-level 

XR(N,I) 

the  quantity  A measured  on  the  side  of  decreasing  z (1*2) 
and  on  the  side  of  increasing  z (1-4),  at  the  local 
development  point. 

YDELTA(N) 

boundary  layer  thickness  at  each  z-station. 

Y (M,N) 

curvilinear  coordinate  normal  to  body  surface. 

2.  Description  of  Subroutines 

The  flow  diagrams  in  Figure  A1  through  A4  describe  the  sequence 
of  events  in  the  computer  code.  Basically,  the  code  is  built  around  a 
driver  program  called  MAIN,  which  sequences  most  of  the  activity  given 
in  the  diagrams.  All  computations  are  performed  through  a series  of 
subroutines  and  functions.  Following  is  a brief  narrative  on  each  of 
these . 

BOUND  (L). 

At  each  x-station  computes  local  edge  conditions  for  each  z-station 
and  local  pressure  gradients  (transformed)  in  the  x and  z directions. 
Linearly  extrapolates  boundary  layer  thickness  from  previous  x-station 
for  each  z.  Forms  a "best  guess"  set  of  profiles  to  start  Iteration 
process  at  the  current  x-station  by  scaling  the  converged  solution  from 
the  previous  x-station  (entry  point  GUESS). 


CONT. 


Computes  boundary  layer  thickness,  maximum  turbulence  Intensity, 


and  wall  shear  stress  for  each  z,  on  each  Iteration.  Solves  continuity 
equation  for  normal  velocity  V. 


DISP(L). 

Computes  displacement  thickness  for  each  z-station  from  the  differ- 
ential equation  given  in  the  text. 


EQTNS. 

Computes  the  4x4  coefficient  matrices  A^  through  A^  given  in 
Equation  (2.32)  through  (2.57).  These  coefficients  are  updated  on  each 
iteration.  This  subroutine  is  called  from  SETUP. 
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FUNCT. 

Forms  a function  used  to  generate  certain  matrix  elements  in 

EQTNS . 

GEOMT(L). 

Loft  data  are  entered  in  this  subroutine  and  all  geometric  func- 
tions are  computed  for  each  surface  point.  The  quantities  are  generated 
at  z+  and  z~  as  appropriate,  and  taking  account  at  the  ends  the  necessary 
values  for  infinite  swept  or  plane-of -symmetry  conditions.  Leading  edge 
and  local  sweep  angles  are  also  computed.  Rotation  matrices  are  gener- 
ated at  each  z-station  for  z+  and  z , for  vector  transfer  in  the  z-direc- 
tion.  Rotation  matrices  are  also  formed  at  each  z for  vector  transfer 
in  the  x-direction.  These  calculations  are  performed  for  each  x-station 
in  succession  and  stored  on  tape.  The  entry  point  GE0M2  then  returns  the 
data  row-wise  for  the  marching  procedure  in  the  main  program. 

INITV. 

Computes  flow  conditions  along  the  swept  leading  edge.  Applies  cri- 
teria to  determine  whether  turbulent  or  laminar  flow  conditions  are  to  be 
applied.  Starting  profiles  for  the  vector  F are  generated  accordingly. 

INPUT. 

Reads  input  conditions  from  cards  and  computes  appropriate  con- 
stants of  the  flow. 

INVRT(A.B.ID) . 

Computes  the  inverse  of  a 4 x 4 matrix.  Inverse  of  A returned  in 
B.  ID-0  returns  singularity  condition.  This  subroutine  is  called  in 

i 1 

SOLVE. 

HULT4(A.B.C.JJ). 

Product  of  Mtrlces  A and  B returned  in  C.  Matrices  are  4 x JJ. 

>• 


Called  In  SETUP  and  SOLVE. 


PRINT(L) . 

Prints  solution  profiles  in  x and  z,  stream  and  normal,  forms. 
Computes  and  prints  wall  shear  stresses,  prints  boundary  layer  thick- 
ness and  displacement  thickness.  Also  prints  iteration  history  of  the 
solution. 

SCALE. 

Used  on  first  x-step  away  from  the  leading  edge  to  scale  the  solu- 
tion back  to  the  attachment  line. 

SEPCON(L) . 

This  subroutine  is  called  when,  for  any  z-station,  the  given  pres- 
sure gradients  cause  separation  to  occur.  The  solution  cycle  then 
passes  through  SEPOON  successively  until  a decreased  adverse  gradient  is 
found  which  allows  the  iteration  cycle  to  be  completed  for  the  given 
x-station.  The  local  adverse  gradient  is  decreased  linearly  on  each 
pass. 

SETUP. 

Sets  up  the  tri-diagonal  matrix,  each  of  whose  elements  is  a 4 x 4 
matrix.  EQTNS  is.  called  to  form  the  A^  through  Ag  matrices  at  each 
mesh  point.  These  are  then  loaded  appropriately  in  the  elements  of  the 
tri-diagonal  system. 

SOLVE. 

Performs  inversion  for  the  tridiagonal.  Called  from  SETUP.  Returns 
a solution  of  the  system  of  equations  for  F on  each  iteration. 

STORE. 

Called  on  each  x-step.  Converged  solution  for  current  x-step  is 
stored.  In  addition  this  solution  is  rotated  to  the  coordinate  system 


1 


i 


for  the  next  x-step,  for  each  z-station.  Entry  point  ST0R2  writes  the 
solution  to  tape. 


ZDERV . 

The  function  ZDERV  performs  z-derivatives  of  any  quantity,  scalar 
or  vector,  one-sided  or  centered,  and  which  may  be  defined  at  some  or 
all  of  the  six  local  developments  points  associated  with  the  point  where 
the  derivative  is  taken. 

i 

i 
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Figure  A2  FLOW  CHART  (continued) 
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